% Iterates on G
%
% Jon Steinsson and Emi Nakamura, July 2006
%*************************************************

function [G,rV,F,F2,Q] = EqSolverGeNCalvoPlus(rPi,menucost,menucost2,...
    RMU,alphaCalvo,beta,theta,a,rp,rad,...
    Proba,ProbNgr,rV,Q,G,SectorWeights,tolV,tolQ,MaxItOnG,yesprint)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

NSectors = size(SectorWeights,1);

rpinc = rp(2)-rp(1);

supDifG = 2*rpinc;
GMove = zeros(radgridsize,1);
ii = 1;
while (supDifG >= 0.99*rpinc && ii <= MaxItOnG)
    if yesprint == 1
        fprintf('######## %5.3f #########\n\n',ii)

%         figure(65)
%         plot(rad,G)
    end
    
    G_old = G;
    
    % Bellman function iteration
    for jj = 1:NSectors
        [rV(jj).rV, F(jj).F F2(jj).F2] = VFIterationGECalvoPlus(rPi,...
            menucost(jj).menucost,menucost2(jj).menucost2,alphaCalvo(jj),...
            rV(jj).rV,RMU,beta,a,rp,rad,Proba(jj).Proba,ProbNgr,...
            G,tolV,yesprint,yesprint);
        if ii == 1
            if jj < NSectors
                rV(jj+1).rV = rV(jj).rV;
            end
        end
    end
    
    % Calculate the stationary distribution
    for jj = 1:NSectors
        Q(jj).Q = StaDistGECalvoPlus(Q(jj).Q,alphaCalvo(jj),Proba(jj).Proba,...
            ProbNgr,F(jj).F,F2(jj).F2,G,a,rp,rad,tolQ,yesprint);
        if ii == 1
            if jj < NSectors
                Q(jj+1).Q = Q(jj).Q;
            end
        end
    end    
    
    % Update G
    [G W_SPold supDifG G_resid GMove] = UpdateGMultiSectorCalvoPlus(Q,F,F2,...
        G,alphaCalvo,a,rp,rad,theta,Proba,ProbNgr,SectorWeights,GMove);

    totWeight = zeros(radgridsize,1);
    for jj = 1:NSectors
        totWeight = totWeight + sum(reshape(SectorWeights(jj)*Q(jj).Q,...
            [agridsize*rpgridsize radgridsize]))';
    end
        
    if yesprint == 1
        [supDifG rpinc]
        [rad G (G-G_old) G_resid totWeight W_SPold]
        fprintf('      rad       G     (G-G_old)    G_resid   totWeight   W_SPold\n')
    end
    
    ii = ii + 1;
end


